package bubbleindex;

import static info.yeppp.Core.Add_V64fS64f_V64f;
import static info.yeppp.Core.Add_V64fV64f_V64f;
import static info.yeppp.Core.Multiply_V64fS64f_V64f;
import static info.yeppp.Core.Multiply_V64fV64f_V64f;
import static info.yeppp.Core.Subtract_IV64fS64f_IV64f;
import static info.yeppp.Math.Cos_V64f_V64f;
import static info.yeppp.Math.Log_V64f_V64f;
import static info.yeppp.Math.Sin_V64f_V64f;
import org.apache.commons.math3.util.FastMath;

public class LombScargle {
      
    static final int FREQ_SIZE = 70;
    static final int Q_SIZE = 18;
    static final int H_SIZE = 19;
    static double[] TestFrequencies = new double[FREQ_SIZE];
    static double[] Q = new double [Q_SIZE];
    static double[] H = new double [H_SIZE];
    
    static final double[] logQi = new double[Q_SIZE];
    static final double[] QiM = new double[Q_SIZE];
    static final double[] cOne = new double[Q_SIZE];
    static final double[] cTwo = new double[Q_SIZE];
    static final double[][] powTempVar = new double[Q_SIZE][H_SIZE];
    
    static final float Omega_float = BubbleIndex.Omega;
    static final double Omega_double = (double) BubbleIndex.Omega;
    static final float M_float = BubbleIndex.M;
    static final double M_double = (double) BubbleIndex.M;
    
    public static void Initialize() {

        double StartingPoint = Omega_double / (2 * 3.14159) - 0.2;
        double Increments = 0.01;
        
        for (int i = 0; i < FREQ_SIZE; i++) {
            TestFrequencies[i] = StartingPoint + i * Increments;
        }
        
        StartingPoint = 0.1;
        Increments = 0.05;
        for (int i = 0; i < Q_SIZE; i++) {
            Q[i] = StartingPoint + i * Increments;
        }
        
        StartingPoint = -1.0 * 0.9;
        Increments = 0.1;
        for (int i = 0; i < H_SIZE; i++) {
            H[i] = StartingPoint + i * Increments;
        }
        
        for (int i = 0; i < Q_SIZE; i++) {        
            logQi[i] = FastMath.log(Q[i]) * Omega_double;
            QiM[i] = FastMath.pow(Q[i], M_double);
            cOne[i] = 1.0 - QiM[i] * FastMath.cos(logQi[i]);
            cTwo[i] = FastMath.sin(logQi[i]);

        }
        
        for (int i = 0; i < Q_SIZE; i++) {
            for (int j = 0; j < H_SIZE; j++) {
                powTempVar[i][j] = FastMath.pow((1.0 - Q[i]), H[j]);
            }
        }   
        

    }

    /**
     * @HQDerivative The method calculates the the largest value of the
     * periodogram given all test values of H and Q.
     * @param TimeValues The array holding the number of days till crash
     * @param Coef The array holding the coefficients of the linear regression
     * @param SIZE The size of the data window. i.e. 256 days
     * @return The maximum value of the periodogram
     */
    
    public static float HQ_Derivative(double[] TimeValues, double[] Coef,
                final int SIZE) {
        
        double C_one, C_two, B_prime, C_prime;
        float Temp, value = 0.0f;
                
        double[] LogTimeValues = new double[SIZE];
        double[] tempFour = new double[SIZE];
        
        Log_V64f_V64f(TimeValues, 0, tempFour, 0, SIZE);
            
        System.arraycopy(tempFour, 0, LogTimeValues, 0, SIZE);
        
        tempFour = new double[SIZE];
        double[] tempFive = new double[SIZE];
        double[] tempSix = new double[SIZE];

            
        Multiply_V64fS64f_V64f(LogTimeValues, 0, Omega_double, tempFour, 0, SIZE);

        Cos_V64f_V64f(tempFour, 0, tempFive, 0, SIZE);
        Sin_V64f_V64f(tempFour, 0, tempSix, 0, SIZE);
        
        for (int i = 0; i < Q_SIZE; i++) {
            
            double Qi_M = QiM[i];
            C_one = cOne[i];
            C_two = cTwo[i];
            double[] tempSeven = new double[SIZE];
            double[] tempEight = new double[SIZE];
            double[] tempNine = new double[SIZE];
            Multiply_V64fS64f_V64f(tempFive, 0, C_one, tempSeven, 0, SIZE);
            Multiply_V64fS64f_V64f(tempSix, 0, C_two, tempEight, 0, SIZE);
            Add_V64fV64f_V64f(tempSeven, 0, tempEight, 0, tempNine, 0, SIZE);
             
            for (int j = 0; j < H_SIZE; j++) {
                
                double tempVar = powTempVar[i][j];

                B_prime = 1.0 * (Coef[1]) * (1.0 - Qi_M) * 1.0 / tempVar;

                //Took out the negative sign for B_Prime

                C_prime = (Coef[2]) * 1.0 / tempVar;
                
                double[] HQDerivativeData = new double[SIZE];
                
                for (int k = 0; k < SIZE; k++) {
                    
                    HQDerivativeData[k] = FastMath.pow(TimeValues[k], 
                            M_double - H[j]) *
                            (B_prime + C_prime * tempNine[k]);
                }
                
                float[] SpectralDensity = new float[FREQ_SIZE];
                
                ComputeLombScargle(LogTimeValues, HQDerivativeData,
                        SpectralDensity, SIZE);
                
                Temp = MaxOneDim(SpectralDensity, FREQ_SIZE);
                if (Temp > value)
                    value = Temp;
            }
        }
        
        return value;
    }
    
    /**
     * @MaxOneDim The method calculates the maximum value in a 1-d array
     * @param Array The array containing the values
     * @param SIZE The size of the array
     * @return The maximum of the array
     */
    public static float MaxOneDim(float[] Array, final int SIZE) {
        float MaxValue = 0.0f;
        for (int i = 0; i < SIZE; i++) {
            if (Array[i] > MaxValue) {
                MaxValue = Array[i];
            }
        }
        return MaxValue;        
    }
    
    /**
     * ComputeLombScargle
     * @param TimeValues
     * @param TimeSeries
     * @param SpectralDensity
     * @param SIZE 
     */
    public static void ComputeLombScargle(double[] TimeValues,
            double[] TimeSeries, float[] SpectralDensity,
            final int SIZE) {
        
        double Omega, SinSum, CosSum, CosResidSum, SinResidSum,
                Tau = 0;
        
        double[] TwoOmegaT = new double[SIZE];
        double[] OmegaTMinusTau = new double[SIZE];
        double[] Residual = new double[SIZE];
        
        //Find mean of Time Series
        double Sum = 0;
        for (int i = 0; i < SIZE; i++) {
            Sum = Sum + TimeSeries[i];
        }
        
        double Mean = -1.0 * Sum * 1.0 / SIZE;
        
        //TimeSeries becomes the residual
        //Subtract_IV64fS64f_IV64f(TimeSeries, 0, Mean, SIZE);
        Add_V64fS64f_V64f(TimeSeries, 0, Mean, Residual, 0, SIZE);
       // for (int i = 0; i < SIZE; i++) {
         //   Residual[i] = TimeSeries[i] - Mean;
        //}
        
        for (int i = 0; i < FREQ_SIZE; i++) {
            
            Omega = 2.0 * 3.14159 * TestFrequencies[i];
            SinSum = 0.0; CosSum = 0.0; 
            
            double tempValue = 2.0 * Omega;
            double[] tempArrayOne = new double[SIZE];
            double[] tempArrayTwo = new double[SIZE];
            
            Multiply_V64fS64f_V64f(TimeValues, 0, tempValue, TwoOmegaT, 0, SIZE);
            Cos_V64f_V64f(TwoOmegaT, 0, tempArrayOne, 0, SIZE);
            Sin_V64f_V64f(TwoOmegaT, 0, tempArrayTwo, 0, SIZE);
            
            for (int j = 0; j < SIZE; j++) {
                //TwoOmegaT[j] = 2.0 * Omega * TimeValues[j];
                //SinSum = SinSum + FastMath.sin(TwoOmegaT[j]);
                //CosSum = CosSum + FastMath.cos(TwoOmegaT[j]);
                SinSum = SinSum + tempArrayTwo[j];
                CosSum = CosSum + tempArrayOne[j];
            }
            
            Tau = -1.0 * FastMath.atan2(SinSum, CosSum) * 1.0 / (2.0 * Omega);
            CosResidSum = 0.0; SinResidSum = 0.0; SinSum = 0.0; CosSum = 0.0;
            
            double[] tempArrayThree = new double[SIZE];
            double[] tempArrayFour = new double[SIZE];
                     
            Add_V64fS64f_V64f(TimeValues, 0, Tau, tempArrayOne, 0, SIZE);
            Multiply_V64fS64f_V64f(tempArrayOne, 0, Omega, OmegaTMinusTau, 0, SIZE);
            
            Cos_V64f_V64f(OmegaTMinusTau, 0, tempArrayOne, 0, SIZE);
            Sin_V64f_V64f(OmegaTMinusTau, 0, tempArrayTwo, 0, SIZE);
            
            Multiply_V64fV64f_V64f(tempArrayOne, 0, Residual, 0, tempArrayThree, 0, SIZE);
            Multiply_V64fV64f_V64f(tempArrayTwo, 0, Residual, 0, tempArrayFour, 0, SIZE);
            
            for (int j = 0; j < SIZE; j++) {
                //OmegaTMinusTau[j] = Omega * (TimeValues[j] - Tau);
                CosResidSum = CosResidSum + tempArrayThree[j];
                //CosResidSum = CosResidSum + Residual[j] * 
                //        FastMath.cos(OmegaTMinusTau[j]);
                SinResidSum = SinResidSum + tempArrayFour[j];
                //SinResidSum = SinResidSum + Residual[j] *
                  //      FastMath.sin(OmegaTMinusTau[j]);
                SinSum = SinSum + FastMath.pow(tempArrayTwo[j], 2.0);
                CosSum = CosSum + FastMath.pow(tempArrayOne[j], 2.0);
                //SinSum = SinSum + FastMath.pow(FastMath.sin(OmegaTMinusTau[j]), 2.0);
                //CosSum = CosSum + FastMath.pow(FastMath.cos(OmegaTMinusTau[j]), 2.0);
            }
            
            CosResidSum = FastMath.pow(CosResidSum, 2.0);
            SinResidSum = FastMath.pow(SinResidSum, 2.0);
            
            SpectralDensity[i] = (float) (CosResidSum * 1.0f / CosSum + 
                    SinResidSum * 1.0f / SinSum);
        }
    }
}